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ABSTRACT 


It  was  the  objective  of  this  investigation  to  study  the  subsonic  flow  over  a  missile 
forebody.  To  this  end  the  NASA  Ames  developed  OVERFLOW,  three  dimensional, 
Navier-Stokes  Code  was  applied  and  detailed  results  were  obtained  for  the  missile 
forebody.  The  flow  conditions  were  set  at  a  Mach  Number  of  0.3  and  Reynolds  Number 
of  two  million  for  each  angle  of  attack  at  0,  2,  6,  10  and  14  degrees.  These  viscous  flow 
solutions  were  compared  with  three  inviscid  flow  solutions,  namely  Slender  Body  Theory, 
von  Karman's  Method  and  the  NASA  Ames  Panel  Code  called  PMARC.  As  expected,  the 
OVERFLOW  and  PMARC  solutions  were  in  good  agreement  at  zero  and  small  angles  of 
attack.  As  the  angle  of  attack  was  increased,  the  OVERFLOW  solution  showed  the 
development  of  the  vortical  flow  separation  on  the  leeward  side  of  the  body. 
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L  INTRODUCTION 


Over  the  years,  a  number  of  methods  have  been  developed  to  provide  the  missile 
aerodynamicist  with  tools  to  predict  the  aerodynamic  characteristics  of  various  missile 
configurations.  They  range  from  purely  empirical  methods  and  simple  slender  body  theory 
to  panel  methods  and  Navier-Stokes  Methods.  A  recent  review  of  these  methods  can  be 
found  in  Reference  1.  This  review  reveals  the  fact,  that  the  prediction  of  strong  viscous 
and  separated  flow  effects,  presents  a  great  challenge. 

Modem  missiles  tend  to  have  small  lifting  surfaces.  Therefore,  the  contributions  of 
the  body  to  the  forces  and  moments,  and  the  understanding  of  the  wing-body  interference 
becomes  very  important.  In  addition,  the  missiles  are  launched  at  relatively  high  angles  of 
attack,  further  exacerbating  the  problem,  due  to  the  vortex  shedding  from  the  missile 
body.  As  a  result,  there  is  a  growing  need  to  compute  the  detailed  flow  separation  and  the 
vortical  flow  features,  which  occur  over  typical  modem  missile  configurations,  using  the 
full  viscous  flow  equations. 

At  NASA-Ames  a  considerable  amount  of  work  was  done  to  study  the  viscous 
flow  surfaces  over  blunt-nose  bodies  of  revolution.  This  work  used  the  NASA-Ames 
developed  OVERFLOW  Code  for  the  solution  of  the  thin-layer  Navier-Stokes  equations. 
References  2  through  5  document  this  effort.  Recently  Reference  6  extended  this  work,  by 
computing  of  the  flow  field  about  a  complete  missile  at  high  angle  of  attack.  The  results  of 
Reference  6  were  compared  to  experiments,  described  in  Reference  7. 

As  pointed  out,  for  example  in  Reference  8,  vortices  are  shed  from  the  forward 
body  and  the  canards.  These  vortices  affect  the  pressure  distribution  and  the  forces  and 
moments  on  the  missile.  Therefore,  in  this  thesis  the  missile  forebody  of  Reference  6  is 
investigated  in  greater  detail. 

Starting  with  the  solution  for  zero  angle  of  attack,  the  angle  of  attack  is  gradually 
increased  to  14  degrees.  This  allows  the  comparison  of  the  Navier-Stokes  solutions  with 
the  simpler  and  faster  methods  for  the  more  benign  flows  at  small  angles  of  attack,  where 
inviscid  methods  can  be  expected  to  give  reasonable  results.  This  comparison  allows 
assessing  the  range  of  validity  of  the  simpler  methods.  Slender  Body  Theory,  von 
Karman's  Method  and  NASA-Ames  panel  code  PMARC  were  chosen  to  compare  with  the 
Navier-Stokes  calculations  from  the  OVERFLOW  Code.  These  methods  are  briefly 
described  in  Chapters  2  through  4.  Chapter  5  presents  the  major  results  of  this  thesis, 
while  Chapter  6  gives  the  conclusions  and  recommendations  for  future  work. 


1 


2 


H.  SLENDER  BODY  METHOD  FOR  APPROXIMATING  FLOW  PAST  BODIES 
OF  REVOLUTION  AT  ZERO  DEGREE  ANGLE  OF  ATTACK 

A.  THEORY 

The  Slender  Body  Method  is  based  on  the  assumptions  of  incompressible,  inviscid 
and  irrotational  flow.  The  governing  equation  is  Laplace's  equation  for  the  potential 
4>(x,r,0)  [Reference  9], 


32<t>  .  d2<t>  .  1  90  .  1  d2(t>  n 

dx2  dr2  r  dr  r2  302 


(2.1) 


Equation  2.1  is  in  cylindrical  coordinates  (x,r,0).  The  body's  angle  of  attack  is 
assumed  to  be  zero,  and  the  body's  longitudinal  axis  is  assumed  to  coincide  the  x  axis. 

Source  superposition  along  the  x  axis,  then  produces  the  well-known  equation  for 
the  disturbance  potential, 


♦  ~H(- 

V  4*Jf 


or 


f(4)4 

+«-*)T 


(2.2) 


(U)  is  the  free  stream  speed.  (4)  is  the  source  location  on  the  body’s  longitudinal  axis.  The 
body  extends  from  0  to  / .  (F )  is  the  change  of  the  body's  cross  sectional  area  along  the 
longitudinal  axis.  It  is  given  by, 


F(x)  =  2nR(x)^^- 
dx 


The  body  radius  is  given  by, 


r  =  R(x) 


(2.3) 


(2.4) 
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Depending  on  the  complexity  of  the  surface  function,  the  integrals  in  Equations  2.5 
and  2.6  can  be  difficult  to  integrate  analytically.  Therefore,  the  MAPLE  Program  was  used 
to  symbolically  integrate  Equations  2.5  and  2.6  and  solve  Equation  2.7. 

B.  COMPUTER  IMPLEMENTATION  OF  THE  METHOD 

MAPLE  is  a  general  purpose,  commercially  available,  computer  code  that  can 
perform  symbolic  mathematics  [Reference  11].  The  code  may  also  be  used  for  numerical 
mathematics.  To  demonstrate  this  program's  usefulness  in  potential  flow  theory,  Appendix 
D  lists  the  implementation  of  Equations  2.5,  2.6  and  2.7  in  Sections  1  through  3 
respectively.  The  MAPLE  Procedure,  given  in  these  Sections,  models  the  flow  over  a 
spindle  body  of  revolution  as  shown  in  Figure  1.  The  spindle  is  at  zero  degree  angle  of 
attack. 


4 


Figure  1.  Spindle  Body  of  Revolution  Surface  Function  and  Shape. 


The  geometry  is  non-dimensional,  so  that  the  length  is  one  and  maximum  diameter 
(tau)  is  16%  of  the  body's  length.  The  maximum  diameter  occurs  at  the  middle  of  the 
body.  Velocity  components  and  Cp  values  for  10  points  on  the  spindle's  surface  are  given 

in  Section  3  of  Appendix  A. 

C.  MAPLE  PROCEDURE  VERIFICATION 

The  parabolic  spindle  shown  in  Figure  1  is  described  by, 

R(x)  =  2tau(x  -  x2)  (2.8) 

Its  first  derivative  becomes, 

^^-  =  2tau{l-2x)  (2.9) 

dx 

The  absolute  magnitude  of  Equation  2.9  will  always  be  less  than  one,  since  the  maximum 
diameter  is  0.16.  This  satisfies  the  slender  body  approximation. 

To  verify  the  spindle  results  in  Section  3  of  Appendix  A,  Reference  10  gives  an 
approximation  of  Equation  2.5  of  the  form, 
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In  this  equation  (F)  is  Equation  2.3,  while  (F  ),  (F  )  and  (F'v)  represent  second,  third 
and  fourth  derivatives  of  the  body's  cross  sectional  area  respectively.  Section  4  of 
Appendix  A  determines  velocity  values  for  Equation  2. 10.  These  are  compared  to  values 
derived  from  Equation  2.5,  using  the  surface  function  for  the  spindle  in  Figure  1.  Figure  2 
shows  the  comparison  between  Equations  2.5  and  2.10  for  48  surface  points.  The 
agreement  is  seen  to  be  excellent.  Both  solutions  exhibit  the  well-known  singular  behavior 
of  slender  body  theory  at  the  nose  and  tail  of  the  body. 


Figure  2.  Comparison  of  Equation  2.5  and  2.10. 
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D.  SUMMARY 


This  singularity  behavior  is  the  result  of  the  integration  of  Equations  2.5  and  2.6  in 
Sections  1  and  2  of  Appendix  A  respectively.  In  both  cases  integration  yields  terms 
containing  the  logarithmic  expression, 

ln(2Vx2+r2  -2x)  (2.11) 

At  the  spindle's  leading  or  trailing  tip  the  radius  approaches  zero.  At  the  tips  the  x 
coordinate  approaches  either  zero  or  one  respectively.  Under  these  conditions  the 
logarithmic  expression  becomes  indeterminate.  As  depicted  in  Figure  2,  the  velocity 

components  do  not  have  a  finite  value  at  the  tips  of  the  body.  Expression  2.11  also 
appears  in  Equation  2.7  (see  Section  3  of  Appendix  A),  indicating  that  the  Cp  at  the  tips 

will  not  be  the  expected  finite  value  of  one.  Another  indication  of  this  indeterminate 

behavior  is  given  in  Sections  5  and  6  of  Appendix  A. 

In  Section  5  as  r  and  x  approach  zero,  the  Cp  value  peaks  at  approximately  0.89 

The  values  decrease  and  do  not  approach  the  expected  Cp  value  of  one.  In  Section  6  as  r 
approach's  zero  and  x  approach's  one,  the  Cp  value  oscillates  about  0.64,  then  suddenly 

decreases  at  the  x  coordinate  of  0.99927.  At  this  point  the  MAPLE  Program  reports  a 
singularity  in  one  of  the  logarithmic  expressions  occurring  in  the  solution  of  Equation  2.7. 

This  implies  that  the  Slender  Body  Method  is  inherently  inaccurate  near  the  leading 
and  trailing  tips  of  a  body  of  revolution.  The  method  can  be  used  to  confirm 
computational  fluid  dynamic  (CFD)  models  in  the  middle  portion  of  the  body,  provided 
that  the  surface  function  is  relatively  simple  and  continuous.  When  surface  functions  can 
not  be  easily  formulated  for  complex  shapes,  the  Panel  Method  might  prove  more  useful. 
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IH.  VON  KARMAN  S  METHOD  FOR  APPROXIMATING  FLOW  PAST 

BODIES  OF  REVOLUTION 


A.  THEORY 


von  Karman's  Method  assumes  a  potential  flow  about  a  body  of  revolution 
[Reference  12].  The  flow  is  considered  incompressible,  inviscid  and  irrotational.  The 
governing  equation  is  Laplace's  equation, 


a>  ay  iay  1  ay 

dz2  dr2  r  dr  r2  362 


(3.1) 


This  represents  the  flow  within  a  body  of  revolution  in  cylindrical  coordinates  (z,r,0) 
[Reference  13].  The  solution  to  this  equation  may  be  formed  by  superposition  of 
elementary  flows.  These  are  represented  by  stream  functions  (\|f).  The  stream  function  for 
uniform  flow  of  speed  (U)  z  direction  is  given  by, 

y  =  ±Ur2  (3.2) 


To  obtain  a  solution  for  flow  over  a  body  of  revolution  at  zero  angle  of  attack  (non-lifting 
flow),  whose  axis  coincides  with  the  z  axis,  it  is  sufficient  to  distribute  sources  along  the 
body  axis.  The  stream  function  of  a  source  located  in  the  origin  of  the  coordinate  system 
is, 


¥  =  - 


mz 


(3.3) 


The  change  in  stream  function  (zA|r)  at  any  point  in  the  flow  from  sources  along  a 
small  interval  (dq)  there  is, 


qM-<;)dq 

ft^+Xz-cf 


(3-4) 
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The  interval  is  measured  along  the  body's  longitudinal  axis.  (#;)  is  the  source  strength  per 

unit  length.  Integration  of  Equation  3.4  along  this  axis  yields  the  induced  stream  function 
at  any  point  in  the  flow.  This  function  is  created  by  sources  on  the  longitudinal  axis.  The 
effect  of  all  the  source  strengths  and  uniform  flow  are  summed  to  give  the  total  stream 
function  at  any  point  in  the  flow.  From  Reference  13  the  summation  can  be  represented  as, 


V(r,*)  =  }w-2+£ 

Z  i 


^/r2+(z-^.+1)2  -J72  +(z-z-)2  q}  (3.5) 


The  (z}+1)  and  (zj)  coordinates  define  the  location  of  the  source  of  constant  strength  {q}) 

on  the  longitudinal  axis.  The  body  of  revolution  is  defined  by  Equation  2.4.  Realizing  that 
the  body  surface  represents  the  stagnation  streamline  \\r  =  0,  equating  the  right  hand  side 
of  Equation  3.5  to  zero  provides  a  means  to  determine  the  unknown  source  strength  q} . 

For  finite  length  bodies,  the  total  sum  of  the  source  strengths  must  be  zero.  This 
additional  condition  must  be  appended  to  the  system  given  by  Equation  3.5.  The  system 
therefore  becomes, 


{const}  =  [dd]{q} 


(3.6) 


The  vector  {const}  is  only  a  function  of  the  known  body  radius.  The  coefficient  matrix 
[dd]  is  a  function  of  the  known  body  geometry  and  source  locations,  while  the  vector  {q} 
is  the  unknown  source  strengths.  The  resulting  system  of  algebraic  equations  may  now  be 
solved  by  common  equation  solves  ,  such  as  Cramer's  Rule,  Gaussian  Elimination,  etc. 

Once  the  source  strengths  are  known  on  the  body’s  surface,  the  flow  velocity 
components  (ur,  uz)  can  be  defined  through  the  first  derivative  of  Equation  3.5.  From 
Reference  13  the  first  derivative  with  respect  to  z  yields  the  radial  component  ur,.  The 
longitudinal  component  uz  is  obtained  from  the  first  derivative  with  respect  to  r.  Using 
Bernoulli's  Equation,  the  non-dimensional  coefficient  of  pressure  (Cp)  is  obtained  as. 


(3.7) 


The  Cp  is  the  ultimate  result  of  this  Method.  Since  the  Cp  is  a  function  of  the  body's 
geometry,  the  value  will  differ  at  various  points  in  the  outer  flow  field  and  on  the  body's 
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surface.  This  condition  makes  the  Cp  parameter  an  excellent  metric  to  make  comparisons 
of  results  from  various  analytic  and  computer  code  techniques. 

B.  COMPUTER  IMPLEMENTATION  OF  THE  METHOD 


Reference  13  lists  a  computer  code  that  implements  von  Karman's  Method.  The 

code  was  called  SPHERE.  A  FORTRAN  listing  of  the  code  is  given  in  Appendix  B.  This 
code  predicts  Cp  values  about  a  sphere  of  radius  (a).  Since  the  flow  field  is  symmetrical, 

the  code  only  calculates  Cp  values  at  surface  points  on  the  upper  half  of  the  sphere.  The 

free  stream  velocity  and  the  radius  of  the  sphere  were  assumed  to  be  one.  SPHERE  used 
Cramer’s  Rule  to  solve  the  system  of  equations. 

Appendix  C  and  D  contain  listings  of  the  two  output  files  created  by  SPHERE. 
The  file  called  SPH.OUT  in  Appendix  C  gives  the  surface  point  location  (r;,  Z;), 
calculated  (Cp)  and  exact  Cp  iCpex)  for  various  numbers  of  surface  points  on  the  sphere. 
The  exact  Cp  values  are  known  to  be  given  by, 


1  4  Va) 


V 


(3-8) 


The  derivation  of  Equation  3.8  may  be  found  in  Reference  13.  Calculated  and  exact  Cp 

values  were  determined  for  nine  surface  points  as  a  check  of  the  code.  Figure  3  shows 
that  the  calculated  Cp  values  plot  identically  on  the  exact  Cp  values  for  the  nine  points. 

However,  when  the  number  of  surface  points  is  increased,  the  calculated  Cp  values 

diverge  from  the  exact.  This  is  expected  because  the  matrix  becomes  increasingly  ill- 
conditioned. 
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Figure  3.  Calculated  and  Exact  Cp  Distributions  for  9  Surface  Points  on  a  Sphere. 
C.  ERROR  NORM  ANALYSIS 

Figures  4  and  5  show  how  the  calculated  Cp  distributions  diverge  from  the  exact 

distribution,  as  the  number  of  surface  points  increases.  This  divergent  behavior  comes 
from  the  ill-conditioning  of  the  coefficient  matrix  [dd]  in  Equation  2.6.  The  growth  of  this 
ill-conditioning  can  be  defined  as, 

{error} = [dd]{q}-  {const}  (3.9) 

error  norm  =  .Jx  {error}2  (3.10) 

Appendix  D  lists  the  ERR.  OUT  file  from  the  SPHERE  Code.  This  listing  gives  the 
results  of  Equation  3.10  for  various  numbers  of  surface  points.  Examination  of  this  data 
and  the  data's  plot  in  Figure  6  clearly  shows  the  error  growth,  as  the  number  of  surface 
points  increases  When  the  number  of  surface  points  exceeds  39,  the  matrix  [dd]  becomes 
so  ill-conditioned  that  the  code  cannot  obtain  the  vector  {q}  in  Equation  3.6. 
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Figure  4.  Calculated  and  Exact  Cp  Distributions  for  33  Surface  Points  on  a  Sphere. 
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Figure  6.  Error  Norm  verses  Number  of  Surface  Points  for  a  Spherical  Body  of 

Revolution. 


D.  SUMMARY 

As  shown  by  von  Karman,  this  method  of  representing  incompressible,  inviscid 
flow  over  a  sphere  can  be  extended  to  the  flow  over  arbitrary  bodies  of  revolution  at  zero 
and  small  angles  of  attack.  However,  the  method  is  limited  to  slender  bodies,  because 
source  or  doublet  distributions  are  used  only  along  the  axis.  In  contrast  to  the  slender 
body  theory  described  in  Chapter  II.,  it  has  the  advantage  of  describing  the  flow  over 
blunt-nosed  bodies. 

A  more  general  method  that  eliminates  the  restriction  to  slender  bodies  is  afforded 
by  surface  distributions  of  sources  and  doublets.  Such  a  method  is  described  in  the  next 
chapter. 
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IV.  PANEL  METHOD  FOR  APPROXIMATING  FLOW  PAST  BODIES  OF 

REVOLUTION 


A.  THEORY 


The  Panel  Method  is  also  based  on  potential  flow  theory  and  generally  parallels  the 
development  of  the  Slender  Body  and  von  Karman  Methods.  Reference  14  gives  a  detailed 
mathematical  description  of  the  Panel  Method.  This  reference  also  describes  how  the 
Method  is  implemented  in  the  computer  code  called  PMARC  (Panel  Method  Ames 
Research  Center).  PMARC  assumes  an  incompressible  and  irrotational  flow.  The  code 
does  allow  a  viscous  boundary  layer  on  the  body,  but  this  option  was  not  used.  All  of  the 
PMARC  results  in  this  thesis  were  based  on  inviscid  flows.  The  method  solves  the  three 
dimensional  form  of  Laplace's  equation, 

V20  =  O  (4.1) 


As  shown  in  Reference  14,  application  of  Green's  Theorem  leads  to  the  following 
integral  equation. 


(4.2) 


((i)  and  ((!„,)  are  unknown  doublet  strengths,  on  the  body's  surface  element  (dS)  and  on 
the  wake  (W)  respectively.  The  doublet  strength  (|xf)  is  zero,  if  the  arbitrary  point  ( P )  of 
interest  is  off  the  surface.  When  the  point  is  on  the  surface,  the  strength  is  2k.  (n)  is  the 
outward  surface  normal  on  the  body’s  surface,  while  (r)  is  the  vector  between  P  and  dS. 
n  and  r  are  known  from  the  body's  geometry.  The  source  strength  (a)  is  given  by, 

(4.3) 


On  the  body's  surface  the  normal  velocity  (Vnorm)  is  usually  zero  or  may  be  prescribed.  The 
free  stream  velocity  vector  (VJ  is  known.  Therefore,  the  source  strength  is  known  in 
Equation  4.2. 
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Equation  4.2  may  be  solved  numerically  by  discretizing  the  body's  surface  into 
panels.  These  panels  represent  areas  of  source  and  doublet  strengths,  that  when  known, 
will  define  the  flow  field  over  the  body's  surface.  This  marks  the  only  difference  between 
the  Panel  Method  and  those  given  in  Chapters  II  and  HI.  Previous  methods  located  source 
strengths  inside  the  body  on  its'  longitudinal  axis.  In  PMARC  the  source  and  doublet 
strengths  are  assumed  to  be  constant  over  a  panel. 

The  discretization  process  results  in  a  system  of  algebraic  equations  of  the  form, 

Ns  .  Ns  /  .  Nw  .  v 

+  +  =  °^j=l,Ns  (4-4) 

*=1  <t=l  1=1 

The  summation  terms  in  Equation  4.4  are  analogous  to  the  integrals  of  Equation  4.2.  ( B ) 
and  (Q  become  respectively. 


k  1 

(4.5) 

’  c>  =~2n 

(4.6) 

The  first  and  second  terms  in  Equation  4.4  represent  the  placement  of  surface  panels  on 
the  body.  The  third  term  accounts  for  the  wake  panels  trailing  downstream  from  the  body. 
In  PMARC  this  equation  is  solved  by  an  iterative  technique  for  very  large  matrices.  The 
solution  occurs  in  a  time  stepping  fashion,  so  that  either  steady  or  unsteady  flows  may  be 
calculated.  The  wake  doublet  strengths  become  a  function  of  the  surface  doublet  strengths 
or, 


(4-7) 

This  means  that  the  PMARC  time  stepping  algorithm  effectively  substitutes  Equation  4.7 
into  Equation  4.4  to  yield  the  surface  doublet  strengths  directly. 

These  doublet  strengths  are  differenced  to  obtain  the  tangential  velocity 
components  on  the  surface  of  each  panel.  The  velocity  components  are  then  used  in, 


16 


(4.8) 


'  8^ 

fn(o-n(t-in 

l  t2J 

Ar  J 

Equation  4.8  defines  the  pressure  coefficient  (Cft)  at  each  panel  (k).  (Vk)  is  the  flow 

velocity  magnitude  at  each  panel.  The  third  term  in  Equation  4.8  accounts  for  pressure 
changes  when  the  flow  is  unsteady.  This  term  is  zero  for  steady  flows,  so  that  Equation 
4.8  becomes  identical  to  Equations  2.7  and  3.7.  Because  the  theory  and  PMARC  allow  the 
effects  of  wakes,  a  problem  arises  as  to  how  to  attach  wake  panels  to  bodies  of  revolution 
to  simulate  the  effect  of  flow  separation. 

B.  MODELING  WAKES  FROM  BODIES  OF  REVOLUTION 

PMARC  allows  wakes  to  be  excluded  from  the  model.  When  this  occurs  PMARC 
will  calculate  Cp  values,  but  can  not  calculate  circulation  about  the  body.  In  this  case  total 

aerodynamic  forces,  such  as  lift,  drag,  etc.  are  not  determined.  Alternatively,  when  wakes 
are  included  in  the  model,  they  must  be  defined  as  either  flexible  or  rigid.  Flexible  wakes 
require  that  the  separation  line  be  specified  in  the  input  file.  Rigid  wakes  require  initial 
wake  panels  be  defined  in  regions  where  the  flow  is  expected  to  separate. 

Reference  14  strongly  implies  that  PM  ARC'S  wakes  should  only  be  modeled  when 
attached  to  traditionally  defined  lifting  surfaces,  that  have  sharp,  linear  trailing  edges.  This 
allows  a  precise  definition  of  a  separation  line  or  initial  wake  panel  layout  along  the  edge. 
The  code  may  then  properly  apply  the  Kutta  condition  at  the  edge.  When  wakes  are 
attached  to  bodies  of  revolution,  without  lifting  surfaces,  the  results  appear  to  be 
questionable,  depending  on  whether  the  body  has  a  tapered,  hemispherical  or  flat  base 
[Reference  15].  A  literature  search  of  References  16  through  19  provided  no  guidance  into 
acceptable  wake  modeling  for  bodies  of  revolution.  Reference  18  claimed  to  have  modeled 
wakes  on  a  body  with  no  lifting  surfaces,  but  provided  no  details  of  the  PMARC  input  file. 

In  view  of  this  situation  the  author  did  not  include  wakes  in  any  of  the  PMARC 
flow  models  in  the  thesis.  All  flow  models  were  considered  to  be  steady.  This  meant  that 
the  third  terms  of  both  Equations  4.4  and  4.8  had  a  value  of  zero.  To  evaluate  PMARC's 
ability  to  compute  fully  attached  flows  the  code  was  first  applied  to  the  sphere  and  the 
spindle,  so  that  the  results  could  be  compared  with  the  Slender  Body  and  von  Karman's 
Methods. 
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c. 


VERIFICATION  OF  PMARC  SPHERICAL  BODY  MODEL 


A  spherical  body  was  constructed  for  the  PMARC  Code.  Appendix  E  lists  the 
input  file  for  the  sphere.  In  the  file  the  name  lists  &BINP8,  &BINP8A  and  &BINP8B 
were  set  to  zero  in  the  file.  These  settings  created  steady  flow.  The  parameters  in  the 
&WAKE1  command  were  also  set  to  zero.  This  eliminated  the  wakes  from  the  model. 
Figure  7  shows  the  surface  paneling  layout. 


Figure  7.  Surface  Panel  Layout  for  a  Spherical  Body. 


The  panels  were  positioned  on  the  sphere  to  yield  nine  Cp  values  that 
corresponded  to  identical  radial  and  longitudinal  locations  of  the  exact  Cp  values  given  in 

Reference  13.  All  PMARC  results  shown  in  this  thesis  were  generated  by  version  12.21,  4 

March  1994,  of  the  code.  Figure  8  shows  excellent  comparison  between  PMARC  and 
exact  Cp  values. 
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Figure  8.  PMARC  and  Exact  Cp  Distributions  for  nine  Surface  Points  on  a  Sphere 

without  wakes. 

D.  VERIFICATION  OF  PMARC  SPINDLE  BODY  MODEL  AT 
ZERO  DEGREE  ANGLE  OF  ATTACK 

This  trend  was  also  confirmed  when  the  spindle  body  was  modeled  in  steady  flow, 
without  wake  panels  and  at  zero  degree  angle  of  attack.  A  typical  PMARC  input  file  for 
the  spindle  is  given  in  Appendix  F.  Figure  9  displays  the  PMARC  spindle  model,  using  the 
surface  function  given  in  Figure  1.  In  Figure  9  a  positive  rotation  about  the  y-axis 

produces  a  nose  up  positive  angle  of  attack,  using  the  right  hand  rule.  The  geometry  was 
non-dimensionalized  similar  to  Figure  1.  Cp  values  were  calculated  and  compared  to 

similar  values  from  the  Slender  Body  Method,  described  in  Chapter  n.  Figure  10  shows 
that  PMARC  results  compare  favorably. 
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Figure  9.  Surface  Panel  Layout  for  a  Spindle  Body. 


Figure  10.  PMARC  and  Slender  Body  Cp  Distributions  for  48  Surface  Points  on  a 
Spindle  at  Zero  Degree  Angle  of  Attack  without  Wakes. 
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E. 


VERIFICATION  OF  PMARC  SPINDLE  BODY  MODEL  AT 
FIVE  DEGREE  ANGLE  OF  ATTACK 


The  PMARC  spindle  model  was  exercised  at  both  zero  and  five  degree  angles  of 
attack.  This  was  done  to  verify  the  model's  physical  behavior  from  a  Cp  point  of  view. 

Figure  11,  12  and  13  compare  Cp  distributions  on  longitudinal  lines  along  the  body's 

suction  or  top  surface,  middle  or  mid  plane  surface  and  pressure  or  bottom  surface 
respectively.  These  longitudinal  lines  can,  alternatively,  be  defined  in  Figure  9  as  lying  in 
the  positive  x,  z  plane  for  the  top  surface,  the  positive  x,  y  plane  for  the  mid  plane  surface 

and  the  positive  y,  negative  z  plane  for  the  bottom  surface.  As  expected  for  a  symmetrical 
spindle,  at  zero  degree  angle  of  attack  the  Cp  distributions  are  the  same  longitudinally, 

regardless  of  surface  orientation. 

When  the  spindle  is  elevated  to  an  attack  angle  of  five  degrees,  suction  or  negative 
pressures  develop  on  the  forward  top  surface  in  Figure  11.  As  expected,  the  flow 
accelerates  over  the  forward  top  surface,  creating  negative  pressures.  On  the  aft  top 
surface  the  flow  experiences  an  adverse  pressure  gradient,  slowing  the  flow  and  creating 
positive  pressure. 


Figure  11.  PMARC  Cp  Distributions  for  48  Surface  Points  on  a  Spindle  at  Zero  and 
Five  Degrees'  Angles  of  Attack  on  the  Top  Surface  without  Wakes. 
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Figure  12.  PMARC  Cp  Distributions  for  48  Surface  Points  on  a  Spindle  at  Zero  and 
Five  Degrees'  Angles  of  Attack  on  the  Mid  Plane  Surface  without  Wakes. 


Figure  13.  PMARC  Cp  Distributions  for  48  Surface  Points  on  a  Spindle  at  Zero  and 
Five  Degrees'  Angles  of  Attack  on  the  Bottom  Surface  without  Wakes. 
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Figure  12  indicates  a  relatively  constant  difference  between  Cp  values  at  attack 

angles  of  zero  and  five  degrees  along  the  spindle’s  mid  plane  surface.  At  zero  degree  a 
fluid  particle  experiences  only  body  fixed,  x  and  y  velocity  components,  as  it  travels  along 
the  body's  mid  plane  surface  [see  Figure  9].  The  particle  travels  circumferentially  and 
longitudinally  along  the  body's  side,  when  the  body  is  rotated  to  five  degrees.  In  this  case 

the  particle's  velocity  vector  contains  all  three  velocity  components  in  the  x,  y  and  z 
directions.  The  total  velocity  is  different  for  either  case,  producing  different  Cp 

distributions. 

As  expected,  in  Figure  13  the  flow  slows  on  the  forward  bottom  surface,  creating 
positive  pressures.  On  the  aft  bottom  surface  the  flow  experiences  a  favorable  pressure 
gradient,  creating  negative  pressures. 

F.  SUMMARY 

The  Panel  Method,  as  implemented  in  the  PMARC  Code,  predicts  Cp  distributions 
for  bodies  of  revolution  that  compare  very  well  with  classical  potential  theories.  The  code 
predicts  the  expected  physical  behavior  of  bodies  at  small  angles  of  attack.  Only  steady, 
incompressible,  inviscid  and  irrotational  conditions  should  be  considered  in  the  flow.  More 
investigation  is  needed  to  understand  PMARC's  predictions  very  near  the  tips  of  conical 
shaped  bodies.  Guide  lines  should  be  developed  for  modeling  wakes  on  bodies  of 
revolution. 

Even  with  these  slight  deficiencies,  PMARC  is  a  suitable  tool  to  confirm  pressure 
fields  about  CFD  models,  as  long  as  viscous  effects  are  small.  Before  explaining  how 
PMARC  may  be  used  to  verify  these  models,  a  brief  discussion  of  the  CFD  computer  code 
used  in  this  study  will  be  given. 
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V.  A  CFD  FLOW  MODEL  OVER  A  TYPICAL  MISSILE  FOREBODY 
A.  A  BRIEF  THEORY  FOR  THE  OVERFLOW  CODE 
1.  Description  of  Basic  Fluid  Dynamic  Equations 

Due  to  the  experience  discussed  in  Chapter  I,  the  Navier-Stokes  flow  solver  called 
OVERFLOW  was  chosen  for  this  study.  OVERFLOW  is  currently  in  development  at  the 
NAS  A- Ames  Research  Center. 

From  Reference  20  the  governing  equations  may  be  written  in  the  three  Cartesian 
x,  y  and  z  coordinates.  In  the  vector  form  the  equation  for  conservation  of  the  flow  mass 
becomes, 

^P+V*pV  =  0  (5.1) 

at 

The  conservation  of  flow  momentum  in  compact  form  is, 

pf  =  V-*s  «.2) 

Three  momentum  equations  are  obtained,  when  Equation  5.2  is  expanded  in  each 
Cartesian  direction.  The  conservation  of  energy  equation  in  the  flow  may  be  written  as, 

p^  =  V.(W7-)+(Js|i  (5.3) 

Dt  dXj 

In  these  five  equations  (p),  (t ),  (V),  (e{),  (k)  and  (T)  are  respectively  density,  time,  the 

velocity  vector,  internal  energy,  thermal  conductivity  and  temperature.  The  fluid  stresses 
(G;j)  are  split  into  pressure  and  viscous  terms  as. 
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(5.4) 


a,=/>5,7+4 


duL+dul 

Kdxj  dx,  j 


The  first  term  in  Equation  5.4  is  the  normal  pressure  (p)  and  (di})  is  the  Kronecker  Delta 


Function.  This  term  represents  the  normal  stress  components  on  the  fluid  element.  In  the 

(  3..  A 


second  term  (p)  is  the  molecular  viscosity  and 


dx; 


Kdxj 


is  the  strain  in  the  fluid.  In  the 


third  term  (X)  is  the  bulk  viscosity  and  is  usually  assumed  equal  to  [  p  |  [Reference 


21].  When  this  assumption  is  made,  the  last  two  terms  in  Equation  5.4  may  be  combined 
to  form  the  viscous  or  shearing  stress  components  on  the  fluid  element. 

In  OVERFLOW  the  molecular  viscosity  in  the  flow  field  is  approximated  by 
Sutherland's  Formula  as, 


p  frY'Vr.+s^ 
4,  [tJ  l t+s) 


(5.5) 


(po)  and  (7^)  are  the  reference  flow  viscosity  and  temperature  respectively.  The  parameter 
(S)  is  assumed  to  be  equal  to  194  degrees  Rankine  for  standard  air.  The  turbulent 
viscosity  is  determined  by  the  Baldwin  and  Lomax  Turbulence  Model  in  the  boundary 
layer.  This  turbulent  viscosity  is  added  to  Equation  5.5.  Other  turbulence  models  are 
available  in  the  code  [Reference  20]. 

Equations'  5.1,  5.2  and  5.3  allow  modeling  of  compressible,  viscous  and  vortical 
flows,  removing  the  basic  assumptions  underlying  the  previous  potential  flow  theory 
methods.  Various  mathematical  operations  are  performed  on  Equations  5.1,  5.2  and  5.3. 
This  effectively  places  them  in  a  finite  difference  form  suitable  for  a  computer  solution. 


2.  Description  of  the  Major  Mathematical  Operations 

Equations'  5.1,  5.2  and  5.3  are  placed  in  conservative  vector  form,  in  order  to 
maintain  the  accuracy  of  the  solution.  The  equations  become, 
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(5-6) 


dQ  dF  dG  dH  n 
dt  dx  dy  dz 

(2  =  [p,p«,pv,pw,ef 

F  =  [pw,  pu2  +  p  -  x„  ,puv -  ,  puw  - ,  (e  +  p)u  -  ux^  -  -  vn^  +  qxf 
G  =  [pv,  pvw -  x^ ,  pv2  +  p - ,  pvw  -  zyi ,  (e  +  p)v - ux^  -  vxw  -  wxyz  +  qy  f 

H  =  [pw,pwM-tXJ,pwv-tyz,pw2  +p-X2Z,(e+p)yv-uxxz -vxyi-wxa  +qzJ 

( u ),  (v)  and  (w)  are  the  velocity  components.  (x,y,  i,j  =  x,y,z )  are  the  viscous  stress 

components.  These  stresses  can  be  formed  from  the  last  two  terms  in  Equation  5.4.  (e )  is 
the  total  energy  and  (q{,  i  =  x,y,z )  are  the  heat  transfer  terms.  Equation  5.6  is  non- 

dimensionalized,  so  that  the  different  flows  may  easily  be  compared,  using  parameters, 
such  as  Mach,  Reynolds  and  Prandtl  Numbers.  Solutions  from  non-dimensionalized 
equations  may  also  be  applied  to  a  flow  over  different  size  bodies  of  the  same  shape.  The 
equations  are  transformed  from  the  physical  domain  into  computational  domain.  This 
transformation  is  done,  so  that  the  equations  may  eventually  be  differenced  on  a  uniformly 
spaced  grid  [Reference  22].  Equation  5.6  is  further  simplified  by  the  thin  layer 
approximation,  which  neglects  the  viscous  terms  over  the  body  in  the  streamwise 
direction..  The  equations  can  now  be  cast  into  finite  difference  form.  This  leads  to  a 
system  of  algebraic  equations,  which  are  solved  implicitly  in  OVERFLOW. 

The  system  in  a  very  simple  form  may  be  written  as, 

[lhs]{AQ}  =  {r/w}  (5.7) 

This  is  called  the  delta  form  [Reference  4],  The  matrix  [Ihs]  and  vector  {r/w}  are 
functions  of  known  flow  quantities,  such  as  density,  pressure,  etc.  The  vector  {A<2}  is  the 
unknown  change  in  the  dependent  variables  between  the  current  and  next  time  step.  For 
steady  flows,  when  the  change  in  these  variables  approaches  zero,  the  solution  is 
considered  to  have  converged  to  the  best  approximation  of  the  original  Equations  5.1 
through  5.3. 
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B.  OVERFLOW  MISSILE  FOREBODY  MODEL  DESCRIPTION 


A  flow  field  about  a  missile  forebody  was  constructed  and  exercised  in 
OVERFLOW  to  obtain  an  acceptable  converged  solution.  Figure  14  shows  the  missile 
forebody's  surface,  but  for  clarity  it  does  not  show  the  flow  grid  extending  above  the 
surface.  The  missile  is  oriented  in  an  X  type  configuration  to  the  flow.  The  axis  system  is 
body  fixed  with  the  same  angle  of  attack  definition  given  in  Figure  9.  All  of  the  dimensions 
of  the  grid  have  been  non-dimensionalized  with  respect  to  the  missile's  maximum  diameter. 
The  maximum  diameter  is  one,  while  the  missile's  length  is  approximately  17.5  diameters. 
The  flow  grid  extends  39  diameters  forward  of  the  nose.  The  grid  extends  over  the  missile 
body's  surface  at  the  aft  end  to  75  diameters  and  for  computational  efficiency  only  one  half 
of  the  flow  field  is  modeled.  The  flow  is  assumed  symmetric  about  the  x,  z  plane.  There  is 
no  grid  in  the  negative  y  direction.  There  is  no  grid  aft  of  the  missile  body’s  rear  base. 

The  grid  over  the  missile  forebody’s  surface  is  composed  of  six  sub-grids.  These 
define  the  in-flow  over  the  missile  forebody's  nose,  the  body's  mid  section,  the  top  and 
bottom  canards,  and  the  out-flow  over  the  body's  rear  section.  The  sixth  grid  is  an  outer 
flow  grid  and  covers  the  mid  section.  These  six  grids  provide  874,125  total  grid  points  to 
describe  the  flow  field  over  the  missile's  forebody. 

All  these  grids  overlap  each  other  to  varying  extents.  In  the  overlapped  regions  the 
chimera  scheme  is  used  to  provide  interpolation  and  grid  blanking  information  [Reference 
23].  This  information  is  used  to  allow  OVERFLOW  to  interpolate  the  flow  variables  from 
one  grid  to  the  next.  The  six  grids  and  the  interpolation  information  are  provided  to  the 
code  in  two  large  files  called  GRID.IN  and  FORT.2  respectively.  The  files  were  from 
Reference  24.  In  addition  to  these  files  OVERFLOW  requires  an  input  file  that  describes 
such  things  as  the  number  of  time  steps  desired,  angle  of  attack,  Mach  and  Reynolds 
numbers,  time  step  size,  boundary  conditions,  turbulence  model,  etc.  Appendix  G  lists  a 
typical  OVERFLOW  input  file  for  the  missile  forebody.  A  description  of  the  various  input 
parameters  can  be  found  in  Reference  20.  All  OVERFLOW  results  shown  in  this  thesis 
were  derived  from  version  1.6ao,  8  February  1994,  of  the  code.  The  forebody  model  was 
initially  exercised  at  a  benign  flow  condition  of  zero  angle  of  attack,  in  order  to  check  the 
flow  physics  and  solution  convergence. 
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c. 


OVERFLOW  MISSILE  FOREBODY  FLOW  SOLUTION 
CONVERGENCE  STUDY 


The  initial  flow  field  was  set  for  a  free  stream  Mach  Number  of  0.3,  a  Reynolds 
Number  of  2.0  million  and  zero  degree  angle  of  attack.  Figures  15  through  20  show  the 
plots  of  the  L,  norm  of  the  vector  { AQ }  in  Equation  5.7  at  each  iteration.  The  L,  norm  is 

the  square  root  of  the  sum  of  squares  of  the  vector's  elements.  As  the  norm  approaches 
zero,  the  change  in  the  dependent  variable  also  approaches  zero,  creating  a  converged 
solution.  The  plots  in  these  figures  indicate  an  acceptable  converged  solution,  after  4000 
time  steps. 

D.  PMARC  MISSILE  FOREBODY  MODEL  DESCRIPTION 

A  PMARC  flow  field  of  the  missile  forebody  was  constructed,  using  4836  surface 
panels.  As  discussed  in  Chapter  IV,  wake  panels  and  boundary  layers  were  not  included. 
The  geometry  of  the  PMARC  missile  forebody  in  Figure  21  was  very  similar  to  the 
OVERFLOW  forebody  in  Figure  14.  Some  slight  differences  existed  between  the  canards 
of  the  two  forebodes.  The  PMARC  model  uses  the  same  body  fixed  coordinate  system 
and  angle  of  attack  definition  in  Figures  9  and  14.  Appendix  H  lists  a  typical  input  file  for 
the  PMARC  missile  forebody.  A  description  of  the  various  input  parameters  can  be  found 

in  Reference  14.  The  PMARC  model  was  exercised  at  various  angles  of  attack  and 
compared  to  Cp  distributions  derived  from  the  OVERFLOW  forebody  model. 

Figure  22  compares  typical  Cp  distributions  from  both  the  OVERFLOW  and 
PMARC  models.  This  figure  indicates  that  most  of  the  Cp  changes  occur  in  the  nose 

region  of  the  missile.  This  region  extends  approximately  three  to  five  calibers  aft  the  nose, 

along  the  x  axis.  Figure  23  shows  the  pressure  distribution  over  the  first  three  calibers  and 
the  missile  body  profile.  Similar  Cp  distribution  comparisons  were  made  in  this  region  at 

various  angles  of  attack,  in  order  to  demonstrate  where  the  two  models  begin  to  differ. 
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Figure  14.  OVERFLOW  Missile  Forebody  Surface. 
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top  surface 


Figure  17.  L,  Norm  for  the  Outer  Flow  Grid. 


Figure  18.  Lj  Norm  for  the  Rear  Section  Grid. 
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Figure  19.  L 2  Norm  for  the  Top  Canard  Grid. 


Figure  20.  Norm  for  the  Bottom  Canard  Grid. 
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E.  COMPUTING  TIMES  FOR  OVERFLOW  AND  PMARC 


Typically  OVERFLOW  calculations  required  approximately  10  to  12  cpu  days  on 
the  Naval  Postgraduate  School,  Cray  Y-MP/EL  computer  for  each  converged  solution.  A 
typical  PMARC  solution  required  1  to  2  cpu  days  on  a  Silicon  Graphics  workstation.  Die 
Cray  mass  storage  disk  space  was  accessed  across  the  network  for  the  PMARC  runs, 
since  local  workstation  space  was  very  limited.  A  PMARC  solution  could  have  been 
obtained  much  faster,  if  more  local  disk  space  were  available. 

F.  OVERFLOW  AND  PMARC  CP  DISTRIBUTION 

COMPARISONS  AT  VARIOUS  ANGLES  OF  ATTACK 

Figures  25  through  38  present  Cp  distributions  for  the  top  ,  mid  plane  and  bottom 

surfaces  of  the  nose  regions  of  the  two  missile  forebody  models.  The  models  were 
exercised  at  0,  2,  6, 10  and  14  degree  angles  of  attack.  The  scale  used  in  these  figures  was 
the  same  as  used  in  Figure  23.  Figures  39  and  53  are  enlargements  of  figures  25  through 
38  in  the  nose  region.  These  figures  show  how  OVERFLOW  and  PMARC  results  diverge 
as  the  angle  of  attack  increases. 

At  zero  and  two  degree  angles  of  attack  in  Figures  24  through  29  and  Figures  39 
through  44,  there  are  very  good  agreement  between  the  two  models.  At  six  degree  angle 
of  attack  in  Figure  30  through  32  and  Figures  45  through  47,  some  differences  begin  to 
appear  along  the  mid  plane  surface.  The  differences  in  the  aft  portion  of  the  bottom 
surface  are  significant  in  Figure  32. 

Figures  33  through  35  and  48  through  50  show  Cp  distributions,  when  the  models 

are  at  10  degree  angle  of  attack.  Both  models  are  in  excellent  agreement  within  the  first 
caliber  along  the  x  axis.  However,  Figures  33  through  35  show  differences  between  both 
models  in  the  aft  portion  of  the  nose  region. 

Figures  36  through  38  and  51  through  53  present  the  results  of  the  two  models  at 
14  degree  angle  of  attack.  These  figures  show  agreement  within  the  first  caliber  along  the 
x  axis,  similar  to  the  10  degree  case.  Again  there  are  differences  between  the 
OVERFLOW  and  PMARC  models  at  14  degree  angle  of  attack  in  the  aft  portion  of  the 
missile's  nose. 
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Figure  21.  PMARC  Missile  Forebody  Surface. 


Figure  22.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Full  Missile 

Forebody. 


Figure  23.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Nose. 
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Figure  24.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Top  Surface. 
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Figure  25.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  26.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  27.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Top  Surface. 


Figure  28.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  29.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Bottom  Surface. 


38 


OVERFLOW 


-x —  PMARC 


5  n 


Canard  Leading  Edge 


x-axis  (non-dimensional) 

Figure  30.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Top  Surface. 


Figure  31.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 
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Figure  32.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  33.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Top  Surface. 
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Figure  34.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  35.  OVERFLOW  and  PMARC  C.  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  36.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Top  Surface. 
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Figure  37.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  38.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  39.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Top  Surface. 


Figure  40.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  41.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Zero  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  42.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Top  Surface. 


Figure  43.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  44.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Two  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  45.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Top  Surface. 
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Figure  46.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  47.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Six  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  48.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Top  Surface. 


Figure  49.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  50.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Ten  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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Figure  51.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Top  Surface. 
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Figure  52.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Mid  Plane  Surface. 


Figure  53.  OVERFLOW  and  PMARC  Cp  Distributions  for  the  Missile  Forebody  at 
Fourteen  Degree  Angle  of  Attack  on  the  Bottom  Surface. 
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G.  SUMMARY 


Missile  forebody  models  were  constructed  for  use  in  both  the  OVERFLOW  and 
PMARC  Codes.  Based  on  Cp  distributions,  comparisons  were  made  between  these  two 

models  for  a  low  Mach  number  at  several  angles  of  attack.  Agreement  was  very  good  at 
very  low  angles  of  attack.  At  six  degree  and  higher  angles  of  attack,  differences  between 
these  two  models  began  to  appear. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  thesis  the  NASA-Ames  developed  OVERFLOW  Navier-Stokes  Code  was 
used  to  obtain  solutions  for  Mach  =  0.3  and  Reynolds  =  2.0e+06  flow  over  a  typical 
canard  configured  missile  forebody.  In  the  absence  of  detailed  experimental  flow  and 
pressure  data,  the  Navier-Stokes  computed  pressure  distributions  were  compared  with 
potential  flow  solutions,  i.e.,  Slender  Body  Theory,  von  Karman's  Method  and  the  NASA- 
Ames  panel  code  PMARC.  As  expected,  the  first  two  methods  were  only  of  limited 
applicability,  but  PMARC  provides  good  agreement  for  zero  and  small  angles  of  attack. 
The  PMARC  and  OVERFLOW  computed  pressure  distributions  started  to  deviate  in  front 
and  close  to  the  canards.  This  indicates  the  presence  of  flow  separation  and  vortical  flow 
development.  Further,  work  is  necessary  to  display  the  flow  features  in  more  detail,  so  that 
the  underlying  flow  physics  can  be  better  understood. 
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APPENDIX  A.  MAPLE  PROCEDURE  TO  CALCULATE  Cp  VALUES  ON  A 

SPINDLE 


> 


>  #  This  MAPLE  Procedure  calculates  the  Cp  values  about  a  spindle  body  of 

>  #  revolution  at  zero  degree  angle  of  attack,  using  potential 

>  #  theory.  The  Cp  values  are  determined  In  Sections  1  through  3. 

> 

>  #  Section  4  determines  the  non-dimensional  longitudinal  velocity  of 

>  #  Equation  2.1 0  and  compares  it  to  Equation  2.5  to  confirm  the 

>  #  procedure. 

> 

>  #  Sections  5  and  6  investigate  the  indeterminate  behavior  of  the  Cp  or 

>  #  Equation  2.7,  as  the  spindle’s  radius  approaches  zero  and  the  x 

>  #  coordinate  either  approaches  zero  or  one  respectively. 

> 


> 

>  #  SECTION  1.  Calculate  non-dimensional  longitudinal  velocity  component, 

>  #  using  Equations  2.8, 2.3  and  2.5. 

> 

>  tau:=  0.16; 


T 


.16 


>  R:= 

>  derivative  R:= 


2.0*tau*  ( (xi)  -  (xiA2) ); 

R.=  .320 1 -.320 

diff(  R,  xi ); 


derivative_R  :=  .320  -  .640  Sj 

>eq_2_3:=  2.0  *  pi  *  R  *  derivative_R; 

eq_2_3  :=  2.0  n  ( .320  §  -  320 12)  (.320  -  .640 

>  !ntegrand_eq_2_5:=  ( eq_2_3  *  (x-xi) )  /  ( ( (x-xi)A2  +  rA2  )A(3/2) ); 

eq  2 

((*-!  f+rf 

>  Eq_2_5:=  ( (1)  /  (4.0*pl) )  *  ( irrt(  lntegrand_eq_2_5,  xi=0.0..1 .0  ) ); 

Eq_2_5  :=  .2500000000  ( -.2048000000  n 

(-21.0jc2-6.x%3  J%2  +  15.00  x-  3.000  +  6.x2  %sj%2+  9.x  i2*  9.x3  -3.0  rz  +  %3j%2-  3.  ^%3  J%l) 
+  . 2048000000  n 

(-6.r2  +  9.3?  +  6.xl%lJ^+^  +  9.xri-3.i2%lJ^+^-6.xl-6.x%lJJ+7  +  %lJ^77)/ 


%1  :=ln (2.JxT+i*  -2.x) 


%2  :=  £  -  2.0  x  +  1.00  +  r* 

%3  :=  ln(2.  J%2+  2.0  -  2.  x) 

> _ 

> 

>  it  SECTION  2.  Calculate  non-dimensional  radial  velocity  component, 

>  it  using  Equation  2.6. 
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>  integrand_eq_2_6:=  ( eq_2_3 )  /  ( ( (x-xi)A2  +  rA2  )A(3/2) ); 

_  4-320  £-. 320 12)  (.320-. 640  £) 

integrand _eq  2_p  :=  2.0 - ~ - — 

>  Eq_2_6:=  ( (r)  /  (4.0*pi) )  *  ( int(  integrand_eq_2_6,  xi=0.0..1 .0 ) ); 

Eq_2_6  :=  .2500000000  r  ( -.2048000000  n  ( -4.00  r2  +  4.0xz  +  7.0xr2-5.0i?-2.x2r2  +  2.x4 

-  6.x  ln(2.  J%1  +  2.0  -  2.  *)^-4.  r4  - 1.0*  +  3.  ^  ln(2 .  J%T  +  2.0  -  2.  *)  J%T)  + 

.2048000000  n  (2.  x4  +  r2  -  2.  x2  r2  +  x2  +  3.  r2  ln(2.  Jj?  +  r2  -  2.  x)  Jx2  +  r2 

-  6.  x  ln(  2.  Jx^+r2  -  2.  x )  r2  Jx^T/2  -  4.  i4  -  3.  x  r2  -  3.  x3 )  /  (  r2  Jx^+r2)  )jn 

%1  :=x2-2.0x  +  1.00  +  r2 


>  #  SECTION  3.  Calculate  the  Cp  values, 

>  #  using  Equation  2.7. 

>  Cp:=  (  -2.0*(Eq_2_5) )  -  ( ((Eq_2_6)**2)  +  ((Eq_2_5)**2) ); 

Cp  :=  -.5000000000  (-.2048000000  x 

( -21.0 x2 - 6. x %3  J%2  +  15.00 x - 3.000  +  6. x2  %3  J%2  +  9. x r2  +  9. x3 - 3.0 r2  +  %3  J%2  - 3. z2 %3 J%l) 
/  y%2  +  .2048000000  it 

(-6. /2  +  9.x3  +  6.x?%l  +9.xr2-3. ^%1  Jx2  +  r2  -6.x2-6.x%l  Jx^  +  r2  +%1  Jx2  +  r2)/ 


J*Ti*)/n  -  .06250000000  r2  (-.2048000000 


V 

w(-4.00/a  +  4.0a2  +  7.0xf2-5.0x3-2.x2r2  +  2.x4-6.x%3r2  J%2 -4.  z4  -  1.0x  +  3.  r2  %3  J%2~) 

t*J%2 

. n(2-*4  +  '2-2.x2z2+x2  +  3.z2%lJx2  +  z2  -6.x%lz2Jx2  +  z2-4.z4-3.xz2-3.x2)\2/ 

+  .2048000000 - 3 - 1  - * - 1  / 

rZJxt  +  r2  } 

n2  -  .06250000000  ( -.2048000000  n 

( -21.0 x2 - 6. x %3  J%2  +  15.00 x - 3.000  +  6. x2 %3  J%2  +  9. x z2  +  9. x3  - 3.0 z2  +  %3  J%2  - 3. z2 %3  J%l) 
jj%2  +  .2048000000  n 

(-6.  z2  +  9.  x3  +  6.  x2 %1  Jx2  +  z2  +  9.  x  z2  -  3.  z2  %1  Jx2  +  z2  -  6.  x2  -  6.  x  %1  Jx2  +  z2  +  %1  Jx2  +  /2)/ 

Jx2  +  z2  )2/n2 

%1  :=  ln(2.  Jx2  +  z2  -  2.  x) 

%2:=x2-2.0x  +  1.00  +  z2 

%3:=ln(2.J%2+2.0-2.x) 
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>  withQinalg): 

Warning:  new  definition  for 
Warning:  new  definition  for 


norm 

trace 


> 

>  points:=  10; 

>  #  pointss  48; 

> 

>  increments  (1.0)  /  ( points  +  1 ); 


points  :=  10 


increment  :=  .09090909091 


> 

>  X(1):=  increment; 


X(l)  :=  .09090909091 


>  r(1):=  2.0  *  tau  *  ( (X(1))  -  (X(1)  A2) ); 

r(l):=  .02644628100 

>  non_dim_ux(1):=  evaif(  subs(  x=X(1),  r=r(1),  Eq_2_5 ) ); 

non_dim_ux(l)  :=  -.08235051778 

>  non_dim_ur(1):=  evalf(  subs(  x=X(1),  r=r(1),  Eq_2_6 ) ); 

non_dim_ur(l)  :=  .2511767847 

>  cp(1):=  evalf(  subs(  x=X(l),  r=r(1),  Cp  ) ); 

cp(l)  :=  .09482965073 

>  for  i  from  2  by  1  to  points  do 

> 

>  X(i):=  X(i-1)  +  increment: 

>  r(i):=  2.0*tau*  ( (X(i))  -  (X(i)A2) ): 

>  non_dim_ux(i):=  evalf(  subs(  x=X(i),  r=r(i),  Eq_2_5 ) ); 

>  non~dim_ur(i):=  evalf(  subs(  x=X(j),  r=r(i),  Eq_2_6 ) ); 

>  cp(i):=  ”  evalf(  subs(  x=X(I),  r=r0),  Cp  7); 

> 

>od: 


>#writeto(cp); 

> 

>  for  i  from  0  by  1  to  points  do 

> 

>  if  (i  =  0)  then 

>  printf(  *  ‘  ):lprint0: 

>  printf(  ‘  tau  =  %4.2f ,  tau  ):lprint(): 

>  printf(  ‘  '  ):lprintO: 

>  printf(‘  x  r  ux  ur  Cp  ‘  ):iprintO: 

>  printf(‘  ‘):lprintO: 

>  else 

>  printf(‘  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  X(i),r(i),non_dim_ux(i),non  dim_ur(i),cp(i)):lprintO 

>  fi 

> 

>od; 
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tau 


0.16 


X 

r 

ux 

ur 

Cp 

0.0909 

0.0264 

-0.0824 

0.2512 

0.0948 

0.1818 

0.0476 

-0.0136 

0.1874 

-0.0082 

0.2727 

0.0635 

0.0244 

0.1300 

-0.0663 

0.3636 

0.0740 

0.0458 

0.0766 

-0.0996 

0.4545 

0.0793 

0.0557 

0.0253 

-0.1150 

0.5455 

0.0793 

0.0557 

-0.0253 

-0.1150 

0.6364 

0.0740 

0.0458 

-0.0766 

-0.0996 

0.7273 

0.0635 

0.0244 

-0.1300 

-0.0663 

0.8182 

0.0476 

-0.0136 

-0.1874 

-0.0082 

0.9091 

0.0264 

-0.0824 

-0.2512 

0.0948 

>  #  SECTION  4.  Calculate  the  approximate  non-dimensional,  longitudinal  velocity  component, 

>#  using  Equations  2.10  and  compare  to  Equation  2.5. 

> 

>Fp1x:=  eq_2_3; 

Fplx  :=  2.0  n  ( .320  %  -  .320  $=2)  ( .320  -  .640  J=) 

>Fp2x:=  diff(  eq_2_3,  xi>; 

Fp2x  :=  2.0  Jt  (.320  -  .640  - 1.2800  n  ( 320 1  -  .320 12) 

>  Fp3x:=  diff(  eq_2_3,  xi$2 ); 


Fp3x  :=  -3.8400  n  (.320  -  .640 1) 

>Fp4x:=  diff(  eq_2_3,  xi$3 ); 

Fp4x  :=  2.4576000  jt 

>  tempi  :=  sqrt(  xl  *  (1 .0-xi) ); 

tempi  :=J|  Jl.O-g 

>  temp2:=  (Fpl x)  *  ( (0.5-xt)  /  (xi*(1 .0-xi)) ) ; 

tempi  -  2.0  * (-320  \ '  -320 12)  (-320  -  .640 1)  (.5  -  £) 

1(10-1) 

>  temp3:=  (Fp2x)  *  ( ln(  2.0*temp1 )  - 1.0 ); 

temp3  :=  (2.O  n (.320  -  .640 §)2  - 1.2800  ji  (.320 1  -  .320 12))  (ln(  2.0ji^)-l.O) 

>  temp4:=  (Fp3x/  2.0)  *  ( 0.5-xI ); 


temp4  :=  -1.920000000  n  (.320  -  .640 1)  (.5  - 1) 

>  temp5:=  (Fp4x/48.0)  *  (1.0  +  (4.0*((0.5-xi)**2))); 

tempS  :=  .05119999999  n  ( 1.0  +  4.0  (.5  -  £)2) 

>  EtL.2_1 0:=  (Fp2x/(2.0*pi))*(ln(r))-((1 .0/(2.0*pi))*(temp2+temp3+temp4+temp5)); 

Eq_2_10 .5000000000  ’(•”'>"«»  t?  - 1.2800, (.320{-320^))  !„(,) 

Jl 

2  Q  « (-320  %  -  .320  g2)  (.320  ■  ,640 £)  (.5  -  £) 

Ki-0-6) 


-.5000000000 


l 

\ 
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+  (2.0 n (.320  -  .640  |)2 -  1.2800 n (.320 §  -  .320 |2)) (ln^.oji  JuT%)  -  l.o) 

- 1.920000000 n (.320  -  .640 £)  (.5  - %)  +  .05119999999  it  ( 1.0  +  4.0  (.5  -  |)2)]/ji 


> 

>poirrts:=  10; 

points  :=  10 

>  #  points:=  48; 

> 

>  increments  (1 .0)  /  ( points  +  1 ); 

increment  :=  .09090909091 

> 

>  X(1):=  increment; 

X(l)  :=  .09090909091 

>  r(1):=  2.0  *  tau  *  (  (X(1))  -  (X(1)A2) ); 

r(l)  :=  .02644628100 

>  non_dim_ux(1):=  evalf(  subs(  x=X(1),  r=r(1),  Eq_2_5 ) ); 

non_dim_ux(l)  :=  -.08235051778 

>  approx_eq_2_5(1):=  evalf(  subs(  xi=X(1),  r=r(1),  Eq_2_10 ) ); 

approx_eq_2_5(l)  :=  -.0815224471 

> 

>  for  i  from  2  by  1  to  points  do 

> 

>  X(i):=  X(i-1)  +  increment: 

>  r0):=  2.0*tau*  ( (X(i))  -  QC(i)A2) ): 

>  non_dim_ux(i):=  evalf(  subs(  x=X(i),  r=r©,  Eq_2_5 ) ); 

>  approx_eq_2_5(i):=  evalf(  subs(  xi=X(i),  r=r(i),  Eq_2_10  ) ); 

> 

>od: 

> 

>  #  writeto(  eq36ux ); 

> 

>  for  i  from  0  by  1  to  points  do 

> 

.  >  if  (i  =  0)  then 

>  printf(  ‘  ‘  ):lprint0: 

>  printf(  ‘  tau  =  %4.2f ,  tau  ):lprint(): 

>  printf(  *  *  ):lprint(): 

>  printf(‘  x  r  Eq_2_10  Eq_2  5  ‘):lprint0: 

>  printf(  •  ‘  ):lprint0: 

>  else 

>  printf(‘  %8.4f  %8.4f  %8.4f  %8.4f,X(i),r(0,approx_eq_2_5(i),non_dim_ux(l)):lprint0; 

>  fi 

> 

>  od; 

tau  =0.16 

x  r  Eq_2_10  Eq_2_5 
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0.0909 

0.0264 

0.1818 

0.0476 

0.2727 

0.0635 

0.3636 

0.0740 

0.4545 

0.0793 

0.5455 

0.0793 

0.6364 

0.0740 

0.7273 

0.0635 

0.8182 

0.0476 

0.9091 

0.0264 

-0.0815 

-0.0824 

-0.0141 

-0.0136 

0.0222 

0.0244 

0.0423 

0.0458 

0.0514 

0.0557 

0.0514 

0.0557 

0.0423 

0.0458 

0.0222 

0.0244 

-0.0141 

-0.0136 

-0.0815 

-0.0824 

> 


> 

>  #  SECTION  5-  Verify  that  as  r  goes  to  zero  and  x  goes  to  zero, 

>  #  Cp  values  will  become  indeterminate. 


> 

>  points:=  50; 


points  :=  50 


> 

>  increments  (1.0)  /  ( 10.0  +  1.0  ); 

increment  :=  .09090909091 


>X(1):= 


>  r(1):= 
>cp(1):= 


Increment; 


X(l) 

2.0*tau*((X(1))-(X(ir2)); 


r(l): 

evalf(  subs(  x=X(1),  r=r(1),  Cp ) ); 


cp(l) 


.09090909091 

.02644628100 

.09482965073 


>  for  i  from  2  by  1  to  points  do 

> 

>  X(j):=  XO-1)  *  increment: 

>  r(i):=  2.0*tau*  ( (X(i))  -  (X(i)*2) ): 

>  cp(i):=  evalf(  subs(  x=X(i),  r=r(i),  Cp ) ); 

> 

>od: 

>  #  writeto(  xO ); 

>  for  i  from  0  by  1  to  points  do 

> 

>  if  (i  =  0)  then 

>  printf(  ‘  ‘  ):lprintO: 

>  prirrtf(  *  tau  =  %4.2f ,  tau  ):lprintO: 

>  printf(  ‘  ‘  ):lprint0: 

>  printf(‘  x  r  Cp 4  ):lprintO: 

>  printf( 4  4  ):lprintO: 

>  else 

>  printf(4  %18.16e  %18.16e  %19.16e  4,X(i),r(i),cp(i)):lprintO; 

>  fi 


> 

>  od; 


tau  =0.16 
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X 


r 


cp 


9 . 090909 0909999990e-02 
8. 2 64 4 62 8 0 99 9 9 997 8e-0 3 
7 . 513 14800 899 99985e-04 
6 . 830 134 554000 00 10e-05 
6 . 2092 13230999 9992e-06 
5 . 64 4 7 39 30 100 000 01e-0 7 
5. 1315811830000005e-08 
4. 665073802999999 le-09 
4.2409761850000004e-10 
3 . 855432 894999 9998e-ll 
3.5049389950000002e-12 
3. 186308 1770000003e-13 
2.8966437970000004e-14 
2.6333125429999999e-15 
2.3939204939999999e-16 
2. 17 629 1358000000 le-17 
1.9784466890000002e-18 
1.7985878989999998e-19 
1 .63 5079 9 0800 0000 le-20 
1.4 8 64362 800 0000 00e-21 
1 . 3513057090000002e-22 
1.22 84597 35000 000 le-2 3 
1.1 1678 1577 000 00 Ole-2 4 
1.0 152559 790000002e-2 5 
9 .22959 98 08999 999 9e-27 
8. 39 054 52 809 9999 9 6e-2 8 
7.6277684370000006e-29 
6 .93433494 3000000 le-30 
6.3039408569999999e-31 
5.7308553250000002e-32 
5.2098684770000003e-33 
4.7362440700000003e-34 
4.3056764270000000e-35 
3.9142512970000000e-36 
3.5584102699999994e-37 
3.2349184269999993e-38 
2.9408349340000003e-39 
2 .67 34 86 3 040 00000 le-40 
2 .430442094 999999 6e-41 
2.2094928140000002e-42 
2. 008629 8 3 10 00000 Oe-4 3 
1.8260271189999999e-44 
1. 660024 654 000000 Oe-4 5 
1.5091 1332200 OOOOOe-4 6 
1. 3719212019999997e-47 
1.24720 10930 000000e-4 8 
1. 1338 19 174 99 9999 8e-4 9 
1. 030744704 9999999e-50 
9. 37 040 64 0900 0000 6e-52 
8. 5 1855 12 8 100 000 14e-5 3 


2 . 6446280999999999e-02 
2 . 6227716680000002e-03 
2 . 4024010460000001e-04 
2 . 1854937759999998e-05 
1 . 9 86 9 35897 000000 le-06 
1 . 8 063 1555699999 99e-07 
1.642 1058939999999e-08 
1 . 4 928236 100000000e-09 
1 . 357 11237 89999998e-10 
1 . 2337385260000000e-ll 
1 . 12 15 8 04 7 8 000 00 Ole-12 
1 . 0196 186 169999998e-13 
9 . 26926015000000 18e-15 
8.4266001380000004e-16 
7 ,6605455800000010e-17 
6.9641323460000009e-18 
6 . 3310294039999996e-19 
5 . 7554812759999995e-20 
5 . 2 3225 570599 99992e-21 
4.7565960960000007e-22 
4 . 3241782680000007e-23 
3 . 9310711519999992e-24 
3.5737010460000005e-25 
3.2488191319999990e-26 
2 . 9534719380000003e-27 
2 . 6849744900000005e-28 
2 .440 8 859 00000000 le-29 
2.2189871820000004e-30 
2.0172610740000002e-31 
1 . 8338737040000003e-32 
1.667 1579130000000e-33 
1.5 15598 102 0000002e-34 
1. 377 8 164 57 000000 le-35 
1 . 2525604 150000002e-36 
1 . 13869 12859999998e-37 
1 . 0351738970000002e-38 
9 . 4106717880000007e-40 
8 .555 15617 1999999 le-41 
7 . 7774147040000000e-42 
7 . 0703770039999999e-43 
6.4276154600000007e-44 
5 . 8432867799999995e-45 
5 . 3120788919999998e-46 
4 . 8291626299999998e-47 
4 . 3901478459999999e-48 
3 . 9910434980000007e-49 
3 . 6282213599999993e-50 
3 . 29 8 383056000000 le-51 
2 . 9985300500000003e-52 
2 . 7259364099999994e-53 


9 . 4 82 96507 3000002 0e-02 
3. 6 1368535200 0000 le-01 
5. 352 1760 8200 0000 8e-01 
6.6798813509999999e-01 
7 . 69335 66 65999 9997e-01 
8.4041487029999995e-01 
8. 8 13381859000000 2e-01 
8.9211484799999996e-01 
8. 727454 5 10000000 le-01 
8. 232300 129999999 8e-01 
7.4356853099999998e-01 
6. 33761003000000 16e-01 
4. 9380743 199999994e-01 
3. 237 078 159 9 9 9999 8e-01 
1.23462 1549999999 9e-01 
-1. 069295499 9 9999 97e-01 
-3. 67 4 67 3 009 9 999 9 9 7e-01 
-6.5815109900000002e-01 
-9. 78980939000000 16e-01 
-1. 32 9956 823 0000 00 le+00 
-1 . 7110787580000002e+00 
-2. 12234 6727 00000 05e+00 
-2. 56 37 6 0760 0000 00 le+00 
-3.0353208079999998e+00 
-3.5370269280000000e+00 
-4 .068879087 0000000e+00 
-4 . 63087728 19999993e+00 
-5.2230215270000002e+00 
-5.8453118119999994e+00 
-6.4977481520000007e+00 
-7 . 180330527 00000 06e+00 
-7 .893058952 000 00 05e+00 
-8. 6359334 170000004e+00 
-9.4089539470000005e+00 
-1. 02 12 12 04 80000 00 le+01 
-1. 1045433100000000e+01 
-1. 1908891739999998e+01 
-1.28024 964 3000000 le+01 
-1.372624 7 190 00 0000e+01 
-1.468014 3960000001e+01 
-1.5664186760000002e+01 
-1. 6678375619999997e+01 
-1.77227 1058999999 8e+01 
-1.8797 191539999996e+01 
-1.9901818540000004e+01 
-2. 103659 1570000002e+01 
-2.2201510679999995e+01 
-2 . 3396575810000005e+01 
-2.4621787019999999e+01 
-2.5877144220000002e+01 


> 


> 

>  #  SECTION  6.  Verify  that  as  r  goes  to  zero  and  x  goes  to  one. 
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Cp  values  will  become  indeterminate. 


>  X(1):=  0.9500; 

X(  1 )  :=  .9500 

>  r(1):=  2.0Mau*((X(1))-(X(ir2)); 

r(l)  :=  .01520000000 

>  cp(1):=  evalf(  subs(  x=X(1),  r=r(1),  Cp ) ); 

cp(l):=  .1752243365 

>  X(2):=  0.999900; 

X(2)  :=  .999900 

>  r(2):=  ,  2.0  *tau*(  (X(1))  -  (X(1)A2) ); 

r(2)  :=  .01520000000 

>  cp(2):=  evalf(  subs(  x=X(1),  r=r(1),  Cp ) ); 

cp(2)  :=  .1752243365 


>  points  :=  29; 


>  increments  o.oooooi ; 


>  for  I  from  3  by  1  to  points  do 


points :=  29 


increment  :=  .1 10'5 


>  X(i):=  XO-1)  +  Increment: 

>  r0):=  2.0*tau*  ( (X(i))  -  (X(i)A2) ): 

>  cp(i):=  evalf(  subs(  x=X(i),  r=r(i),  Cp ) ); 

> 

>  od: 

Error,  (in  evalf/ln)  singularity  encountered 

>#writeto(x1 ); 

>  for  I  from  0  by  1  to  points  do 

> 

>  if  (i  =  0)  then 

>  printf(  *  '  ):lprint(): 

>  printf(  ‘  tau  =  %4.2f,  tau  ):lprint(): 

>  printf(  ‘  ‘  ):lprint0: 

>  printf(‘  x  r  Cp  ‘  ):lprintO: 

>  printf(  •  •  ):lprint0: 

>  else 

>  printf(‘%18.16f  %18.16f  %19.16f  \  X(i),r(i),cp(i)):lprirrtO; 

>  fi 


tau  =0.16 


0.9499999999999998  0.0152000000000000  0.1752243365000000 
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0.9998999999999999 

0.9999009999999999 

0.9999020000000001 

0.9999029999999999 

0.9999040000000001 

0.9999050000000002 

0.9999060000000001 

0.9999069999999999 

0.9999079999999999 

0.9999090000000000 

0.9999100000000000 

0.9999110000000002 

0*9999120000000001 

0.9999129999999999 

0.9999140000000001 

0.9999150000000002 

0.9999160000000001 

0.9999169999999999 

0.9999179999999999 

0.9999189999999999 

0.9999200000000000 

0.9999210000000001 

0.9999219999999999 

0.9999230000000002 

0.9999240000000001 

0.9999250000000000 

0.9999260000000000 

0.9999269999999999 


>  #  file  s  app_a.txt 

> 

> 

> 


0.0152000000000000 

0.0000316768640000 

0.0000313569280000 

0.0000310369920000 

0.0000307170560000 

0.0000303971200000 

0.0000300771840000 

0.0000297572480000 

0.0000294372800000 

0.0000291173440000 

0.0000287974080000 

0.0000284774720000 

0.0000281575360000 

0.0000278375680000 

0.0000275176320000 

0.0000271976960000 

0.0000268777280000 

0.0000265577920000 

0.0000262378560000 

0.0000259178880000 

0.0000255979520000 

0.0000252780160000 

0.0000249580480000 

0.0000246381120000 

0.0000243181440000 

0.0000239982080000 

0.0000236782400000 

0.0000233583040000 


0.1752243365000000 

0.6476606865000000 

0.6553688371999999 

0.6422468258999999 

0.6393732819000001 

0.6585344844000000 

0.6561241774000001 

0.6536068717000000 

0.6398706602999999 

0.6489621221999999 

0.6580446579000001 

0.6431658841000000 

0.6400647958000000 

0.6507920520000000 

0.6343043787000000 

0.6448765211000000 

0.6425597706000000 

0.6392634326000001 

0.6357834041000000 

0.6482360795000001 

0.6293168993000000 

0.6416774725000000 

0.6547907977000000 

0.6354180657000001 

0.6494795412000000 

0.6623045039000000 

0.6253918384000000 

5.8453118119999994 
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APPENDIX  B.  FORTRAN  LISTING  OF  THE  SPHERE  CODE 
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0  Calculate  tha  surface  point*,  vhera,  17  continue 

e  (at  ie  tha  aphara'a  radius,  do  19,  j  •  1,  n 

o  (ii)  la  tha  horisontal  axle  of  tha  ephere,  dd(n,j)  ■  a(J) 

c  jri)  la  tha  perpendicular  diatanea  freai  (si).  19  continue 
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Jui  20 199508-39 
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S.22498S647159489E-O01 
5. 87785252292473  IE-001 
6 . 49448048330163 7E-00 1 
7 . 07 10678 11865476E-00 1 
7 .604059656000309E-00 1 
8.090169943749475E-001 
8.526401643540922E-001 
8 . 9 1006524186367 9E— 001 
9 -238795 325 112 867E-001 
9.510565162951535E-O01 
9 . 723699203976766E-00 1 
9 . 876883405951378E-001 
9.969173337331280E-001 

CP. 

8.805393639067629E-001 
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In  subroutine  CRAMER, 
denou«  O.OOOOOOOOOOOOOOOE+OOO 
causes  <x(fc)>  in  subroutine  to  blow  up, 
CAN  NOT  obtain  a  (q)  vector  at, 

41,  surface  points  ftl 
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APPENDIX  D.  ERR. OUT  FILE  FROM  THE  SPHERE  CODE 
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This  is  the  data  output  file  called  ERROR. OUT 
from  the  program  SPHERE. FOR. 
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7.049768158406932E- 
1.2 04559 7 1973 1763E- 
2 . 9 89 18729209 14 45E- 
2 . 420 6 4 965 46 1035 3E- 
1. 827890796692 02 3E- 
2 . 427334109556196E- 
4. 8 662 0 8540 02 342 7E- 
2 . 49498 16993 0452 5E- 
3 . 456497292250762E- 
2 . 347015589325810B- 
1.638605492134842E- 
1.01482465089 18  06E-> 
6. 80 4 778254335 41 9E- 
6.117182181739560 
2 . 74 08 4 13024 055 3 IE-1 
1. 655699 71275 489 0E-J 
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